Predicting rare events in chemical reactions: application to skin cell proliferation 
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In a well-stirred system undergoing chemical reactions, fluctuations in the reaction propensities 
are approximately captured by the corresponding chemical Langevin equation. Within this context, 
we discuss in this work how the Kramers escape theory can be used to predict rare events in chemical 
reactions. As an example, we apply our approach to a recently proposed model on cell proliferation 
with relevance to skin cancer [P. B. Warren, Phys. Rev. E 80, 030903 (2009)]. In particular, we 
provide an analytical explanation for the form of the exponential exponent observed in the onset 
rate of uncontrolled cell proliferation. 
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I. INTRODUCTION 

Noise is ubiquitous in systems undergoing chemical re- 
actions. In a well-stirred system, the source of noise 
comes from the probabilistic nature of the reactions, 
and can be analyzed try employing the Chemical Master 
Equation (CME) [1 %. With the help of the Kramers- 
Moyal expansion, a Chemical Langevin Equation (CLE) 
can be formulated to approximate the CME [2 4]. When 
the number of molecules in the system is small, the limi- 
tations of the approximation have been explored in @, [g . 
On the other hand, when the numbers of molecules of 
the different chemical species in the system are greater 
than certain thresholds, the CLE constitutes a reason- 
able approximation to the CME [4J. One big advantage 
of the CLE is the well developed analytical tools avail- 
able. For instance, thermally activated escape theory 
(see, e.g., 0, Q), such as Kramers escape theory, serves 
as a natural platform for the studies of extinction rate of 
chemical species [9l4l3l| , and transition rates between two 
metastable states of the system concerned [L4J. This is 
the approach adopted in this work. Besides being of gen- 
eral interest to chemical systems, the method discussed 
here is also relevant to cellular processes. One interesting 
example is the recent proposal that metastability in skin 
cell proliferation constitutes a component in the patho- 
genesis of cancer [15J. In particular, the author in [l5[ 
observed numerically that the rate for the onset of uncon- 
trolled cell proliferation has an exponential component 
that scales in a specific manner with the model param- 
eters. As an illustration, we shall demonstrate how the 
form of the exponent observed can be explained analyti- 
cally within the context of Kramers escape theory. 



II. A SIMPLE EXAMPLE 

We will first start by considering a simple example to 
set up the formalism. Consider the following set of chem- 



FIG. 1: The potential energy landscape corresponding to the 
model described in Eq. fT}. There exist two fixed points: a 
stable one at n A — M and an unstable one at ua — N. If ua 
gets beyond N, it will diverge to infinity. The phase portrait 
of the model is depicted under the plot. 




ical reactions: 



.4 



A + A 



A + A 



A 



(1) 



where A depends on the number of A molecules, n a , in 
the following manner: 



X(n A ) = 



my , 

2Tn A 



if n A < N 
otherwise , 



(2) 



where M and T are constant. In a deterministic system, 
the above scenario is governed by the following equation: 



n A = [\{n A ) -Tn A ]n A 



(3) 
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In other words, if < n A (t = 0) < N, then n A (t — > 
oo) = M. On the other hand, if n A (t = 0) > N, then 
n A (t — > oo) diverges (c.f. the phase portrait under the 
plot in Fig. [J). 



This deterministic picture is of course incomplete due 
to the neglect of the intrinsic fluctuations from the re- 
action propensities. Such fluctuations are approximately 
captured by the following CLE j4j: 

fiA = [A(ru) - IX4] n A + ^X(n A )n A w + ^Tn 2 A w' , (4) 



where - 



are Gaussian noises with zero means and unit 



standard deviations. Since w and w' are uncorrelated, we 
have 



(5) 



n>A = 7/ + v 2 7f" 
where we have introduced the following functions: 



7 = 



/ 



[A(tia) + Tn A ]n A 

2 
2[A(n^)-rn A ] 
A(n j4 ) -t-TnA 



(6) 

(7) 



Note that the form of Eq. ([5J corresponds to a Langevin 
equation describing a particle in a potential well under 
thermal perturbations in the non-inertia regime. Specifi- 
cally, n A can be treated as the coordinate of the particle, 
/ as the force exerted on the particle due to an underly- 
ing potential, and 7 as the position-dependent damping 
coefficient of the system. 

In the region n A < M, the force is 



f(n A ) = 



2(M-n A ) 

M + n A 



(8) 



The corresponding potential energy can thus be deter- 
mined as 



U(n A ) 



dxf(x) 



(9) 



m 



= 2(ru-M)-4Mln^tM. (10) 

Similarly, in the region n A > M, we have f(n A ) = „- and 
U(n A ) = —2n A /3 + constant. The shape of the potential 
energy is depicted in Fig. Q] 

We are primarily concerned with rare escape events 
and we thus assume that N— M ^> 1. In this scenario, if 
the initial state of the system is such that n A (t = 0) = M, 
then the waiting time, r, for the system to move out of 
the potential well, i.e., the waiting time for n A to attain 
the value N is (c.f. Sect. 7.2 in Q): 



2tt 



U"(M) 



exp 



2(N-M) -AM In 



N + M 



2nM exp 



2(N-M) -4Mln 



2M 

N + M 
2M 



(11) 
(12) 



Note that tia will, with probability one, either go to zero 
or diverge [8| , and so the knowledge of the escape rate is 
particularly important. 



FIG. 2: (Color online) The flow lines of the cell proliferation 
model described in Eqs (|13p and (|14p . The two fixed points of 
the models are depicted by the gray circle and square (c.f. Eqs 
(f22)l and ((23j)). The red (gray) wriggly line depicts schemat- 
ically a possible escape trajectory. The numerical values for 
the model parameters are a — 2, fi — 10, r — 0.08, no = 200, 
T = 0.045, p = 0.22 and Pl = 0.26 [TJ]. 
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III. SKIN CELL PROLIFERATION 

We now move onto discussing a model for skin cell 
proliferation. The model we study is based on the single 
progenitor cell model introduced in [lTJ, [l7| • This model 
was then generalized in [TJj to account for the home- 
ostasis of the system. Furthermore, the author in [15[ 
suggests that the escape of the system from the home- 
ostatic basin due to rare stochastic fluctuations plays a 
role in uncontrolled cell proliferation. This is of impor- 
tant relevance to the study of skin cancer. Specifically, 
there are two basal layer cell types in this model: pro- 
genitor cells A and postmitotic cells B. These two types 
of cells proliferate according to following scheme: 



A^hA + A , A^A + B 
A A B + B ,BA|. 



(13) 
(14) 



The first three processes represent the different progeni- 
tor cell division pathways, and the fourth represents post- 
mitotic cells leaving the basal layer. In the model above, 
r is a constant and Ai are defined as follow: 



where 



Ax = A(n)r(l - q(p)) 


(15) 


A 2 = A(n)(l-2r) 


(16) 


A 3 = A(n)r(l + q(p)) 


(17) 


n A 
= n A + n B , p= — 
n 


(18) 



aw = a (2a 

v n 



An = 



r(i-po) 

Po 



and 



q(p) = tanh 



10po(l-po)(p-p )(pi - p) 
p(l-p)(pi ~ Po) 



(19) 



(20) 



Note that the constants n , p Q and p\ represent the ini- 
tial number of cells, the fraction of progenitor cells at the 
stable fixed point (marked by the gray circle in Fig. [5]) , 
and the fraction of progenitor cells at the unstable fixed 
point (marked by the gray square in Fig. [5]), respectively. 
An experimentally motivated set of parameters for this 
model is shown in the caption of Fig. [2J We refer the 
readers to |15l - fl7| for more detailed physiological inter- 
pretations of the different processes. Here, we will only 
note that a divergence of the total cell density, ha + Ub , 
signifies the onset of uncontrolled cell proliferation. In 
[15j . the author employs the standard Gillespie kinetic 
Monte Carlo algorithm [lj| to perform stochastic simu- 
lations of the model, and finds that the escape rate has an 
exponential component of the form exp(— 4.6 x noAp 2 ) 
where Ap = p\ — po denotes the difference in the frac- 
tions of progenitor cells at the saddle point and at the 
fixed point. We will now demonstrate how the Kramers 
escape theory accounts for the exponent observed. 

We shall first look at the deterministic case, where the 
chemical reaction scheme in Eqs (fT5|) and (fl4l lead to 
the following set of ordinary differential equations: 



iiA = —2\rqriA , hs = A(l + 2rq)n J 



Trii 



(21) 



Setting the L.H.S. in the equations above to zero, we find 
two nontrivial fixed points: 



X = (n p , (1 - p )n ) 
Y = 



(22) 




-Wi 
r(i-pi) 



AoMl-ft) ) . (23) 



These fixed points are denoted by a gray circle and a gray 
square respectively in Fig. [5] along with flow lines. 
The corresponding CLE for this system is Q: 



ha = (Ai - X 3 )ua + ^/\in A wi - ^X 3 n A w 3 (24) 

h B = (A 2 + 2X 3 )n A - Tn B (25) 

+ \J\ 2 nAW2 + 2y / \ 3 nAW 3 - \/YnBW± , (26) 

where the w% are again Gaussian noises with zero means 
and unit standard deviations. As among the four inde- 
pendent Gaussian noise terms, only w 3 are common in 
both equations, we can thus simplify the above equation 
to the followings: 



f>B 



(Ai - X 3 )n A + \/(Ai + \ 3 )n A wi 

(A 2 + 2X 3 )n A - Tn B 

+ \/(A2 + 4A 3 )r7, j4 -(- Tn B crwi + y 1 - a 



(27) 
(28) 



c w 2 



FIG. 3: The magnitude of the correlation a as a function of 
riA and ub around the two fixed points indicated again by the 
square and the circle. The broken line denotes the escape path 
corresponding to the analytical calculation in Eq. (|36|) . and 
the solid line denotes the escape path obtained numerically 
(c.f. the discussion before Eq. Ij39| l1. which corresponds to the 
result in Eq. (f3u]l . 




where 



2\ 3 n A 



yVi(Ai + A 3 )[(A 2 + 4A 3 )n A + Tn B ] 



(29) 



corresponds to the correlation between the two fluctua- 
tion processes. 

Note that in dimensions higher than one, one cannot 
in general represent the force fields as the gradients of 
a potential, i.e., the force is not conservative. Although 
a potential energy cannot be constructed here, it is still 
possible to obtain a scalar function that serves to deter- 
mine the exponent in the Arrhenius term associated to 
the escape process [lj|. This can be achieved by solv- 
ing a second-order boundary value problem, and usually 
can only be done numerically. Here, we will avoid this 
numerical challenge and aim to proceed analytically by 
making a series of approximations to the above CLE. 

As aforementioned, the Gaussian noises associated to 
the two coordinates are correlated. In Fig. |31 we show 
the magnitude of a around the two fixed points, which is 
bounded above by 0.31. The first approximation is that 
we will set a to zero, i.e., we assume that the perturba- 
tions acting on ua and ns are uncorrelatcd. With this 
simplification, Eqs ([24| and (|25j) can be written as: 



2j b wb , 
(30) 




, n B = JbJb - 

A(l + 2r + Arq)n A + Tn B 

2 

2A(1 + 2rq)n A - 2Tn B 

A(l + 2r + Arq)n A + Tn B 



(31) 

.(32) 



FIG. 4: The force vector fields, (/a,/s), according to Eqs 
(|32[) . The magnitudes of the vectors are scaled up uniformly 
for visual clarity. 



3.1 x n Ap 2 



(37) 
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Fig. [4] show the force vectors (/a>/b) around the two 
fixed points, which suggests that \fs\ "C |/a| m the re- 
gion connecting X and Y. In other words, it is much 
easier for the particle to diffuse vertically than to dif- 
fuse horizontally. We will therefore ignore the second 
dimension and consider purely the first coordinate. This 
constitutes our second approximation and effectively col- 
lapses the problem into a one-dimensional problem. As 
a result, we can calculate the corresponding potential by 
simply integrating over /a'- 



U= 2q 



n B 



d.r 



(33) 



where we will take Ub — Xb (c.f. Eqs (|2"0f and 
If the initial cell densities are in the metastable region 
around X, the rate, R, at which uncontrolled cell prolif- 
eration occurs will be of the form R ex exp(— AU) [lj|, 

where 



AU = 



Y A 



■2q 



X A 



X + n B 



dx 



(34) 



Note that our consideration effectively amounts to calcu- 
lating the first passage time of a particle constrained to 
diffuse along the horizontal path depicted by the broken 
line in Fig. [3] 

Since throughout the range of the integration, the ar- 
gument in q is bounded above by 0.12, q is well approxi- 
mated by q where (c.f. Eq. ([20)1 ) 



q{p) 



10/Oq (1 - Po)(p - Po){pi - P) 



p(l- p)(p! 



Po) 



(35) 



This simplification allows us to perform the integration 
in Eq. ([M]) analytically, and we find that for small Ap, 



AC/ 



5p (3-2p )< 
3(1 -Po) 



n Ap 2 +O(Ap 3 ) (36) 



where the second approximated equality comes from sub- 
stituting in the numerical values of the parameters shown 
in the caption of Fig. [2] Recall that the exponent is found 
numerically to be (4.6 x no Ap 2 ) [la]. We have therefore 
recovered the scaling of the exponent with respect to Ap. 
On the other hand, the prefactor we obtained is about 
two thirds of that observed from simulations. 

We will now try to incorporate the second dimension 
and the correlation in the fluctuations into the picture. 
Our strategy is to find a path that better represents the 
escape route. In the weak noise limit, such an optimal 
escape path encapsulates the information on the asymp- 
totic behavior of the escape process, and in principle, 
can be obtained by solving a set Hamiltonian equations 
with the appropriate end points. @, [20l - [22J | . We find 
the application of the numerical procedure to the prob- 
lem concerned challenging as the corresponding set of 
Hamiltonian equations are higher sensitive to the initial 
conditions chosen. Hence, we will instead make a crude 
estimate on the escape path that connects the metastable 
state to the saddle point [24| . Specifically, as we have ar- 
gued that the particle diffuses more easily along the verti- 
cal direction, we reason that as the particle goes upward 
in the nA direction, it should stay at the bottom of the 
valley with respect to the force fields along the ha- We 
therefore start at the metastable point, X, and find the 
tia that minimizes \/a\ as we move up in the hb dimen- 
sion. We find that such a path corresponds to a slanted 
as depicted by the solid line in Fig. [3[ When the path 
reaches Yb in the ns coordinate, we simply connects it 
horizontally with the saddle point Y (c.f. Fig. [3]). We 
denote this escape path by z(s), where < s < L cor- 
responds to the parametrization of the curve such that 
z(0) = X, z(i) = Y and |z'(s)| = 1 for all s. Now, we 
collapse again the problem into one dimension by con- 
sidering only fluctuation processes along this path. The 
time evolution of the particle along the path can be ex- 
pressed as: 

8 = U7a/a + V7b/b (38) 

1 lo 

+ [2(u 2 'fA + v 2 j B + 2auv^/"fAlB)} w , 



where the 7a/b and /a/b are as expressed in Eqs ([31]) 
and (l32j) . a is defined in Eq. ([29]), and (u, v) denotes the 
unit tangent of the curve z at the point (riAi^s)- As 
a result, the exponent in the rate describing the escape 
process from s = to s = L is: 



AU = - 



WJa/a + VlBJL 



o w 7A + v 2 jb + 2auVy/%ATB 



d.s- 



The numerical value is found to be 



AU = 2.05 = 6.4 x n Ap 2 , 



(39) 



(40) 



which is greater than the simulation results by about 
39%. The discrepancy here is likely an outcome of 



our crude way of estimating the escape path. To im- 
prove upon this result, more sophisticated numerical ap- 
proaches would be required, which is beyond the scope 
of this work. 



reactions due to stochastic fluctuations. As an applica- 
tion, we have considered a model on cell proliferation and 
explained analytically the observed rate for the onset of 
uncontrolled cell growth. 



IV. CONCLUSION 

In summary, we have discussed how the Kramers es- 
cape theory can be used to predict rare events in chemical 
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